Europhysics Letters 



PREPRINT 



Attraction and unbinding of like-charged rods (*) 

A. Naji1(**), a. Arnold^, C. Holm^ and R. R. Netz^ (**) 

^ Sektion Physik, LMU, Theresienstr. 37, D-80333 Miinchen, Germany. 
^ MPI fiir Polymerforschung, Ackermannweg 10, D-55128 Mainz, Germany. 



PACS. 87.15.-v - Biomolecules: structure and physical properties. 
PACS. 61.20. Ja - Computer simulation of liquid structure. 

PACS. 87.15.Nn - Properties of solutions; aggregation and crystallization of macromolecules. 

Abstract. - We investigate the effective interaction between two like-charged rods in the 
regime of large coupling parameters using both Molecular Dynamics simulation techniques and 
the recently introduced strong-coupling theory. We obtain attraction between the rods for 
elevated Manning parameters accompanied by an equilibrium surface-to-surface separation of 
the order of the Gouy-Chapman length. A contiimous unbinding between the rods is predicted 
at a threshold Manning parameter = 2/3. 

Introduction. - In recent years, both experiments [1] and numerical simulations [2-10] 
showed that like-charged macroions can attract each other via effective forces of electrostatic 
origin; one prominent example is the formation of dense packages of DNA molecules (DNA 
condensates) [1]. These observations indicate like-charge attraction in the presence of multi- 
valent counterions, at low temperatures or for large surface charge density on macroions, i.e. 
when electrostatic correlations between charged particles become increasingly important. The 
standard mean-field Poisson-Boltzmann (PB) theory predicts repulsion between like-charged 
objects [11]. In contrast, incorporation of electrostatic correlations generates attractive inter- 
actions in agreement with a variety of experimental and numerical results [12-22,27]. The 
importance of electrostatic correlations can be quantified by means of the coupling parameter, 
S = which depends on the charge valency of counterions, q, the surface charge 

density of macroions, (Jg, and the Bjerrum length, £b — e^ / {Aireeok-BT) (associated with a 
medium of dielectric constant e and at temperature T). The PB theory is valid in the limit of 
vanishingly small coupling parameter S ^ [21-23], while non- mean-field features emerge in 
the opposite limit of large coupling, S 3> 1, and are typically accompanied by the formation 
of strongly-correlated counterion layers. For charged curved surfaces, one has to consider the 
entropy-driven counterion-condensation process as well. For rod-like macroions, as consid- 
ered in this paper, counterion condensation is controlled by the so-called Manning parameter, 
£, = q^BT [24], where t stands for the single-rod linear charge density (in units of the ele- 
mentary charge e). For very small Manning parameter, ^, counterions de-condense from the 
rods leading to a bare electrostatic repulsion between them [25]. For sufficiently large ^, on 
the other hand, a certain fraction of counterions remains bound to the rods and attraction 
becomes possible for moderate to large couplings. We are interested in the regime of Manning 
parameters, ^, and coupling parameters, S, where an effective rod-rod attraction is present. 

Several mechanisms for the attraction between like-charged rods have been considered 
in the literature, including covalence-like binding [26], Gaussian-fluctuation correlations [14], 
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and structural correlations [18-20]. The two latter approaches yield attraction based on cor- 
related fluctuations of condensed counterions on opposite rods and short-ranged correlations 
due to the ground-state configuration of the system, respectively. Conflicting predictions for 
the threshold value of Manning parameter, above which attraction between two rods is pos- 
sible, have been obtained: The analysis of Ray and Manning [26], based on the standard 
counterion-condensation model [24], leads to attraction for ^ > 1/2. In contrast, an attraction 
regime of ^ > 2 was proposed by the counterion-condensation theory of Arenzon et al. [20] . 
Recent numerical simulations [3,10], on the other hand, reveal attraction already for Manning 
parameters of the order of ^ « 1 and for moderate coupling parameters, though did not specif- 
ically consider the threshold value of ^. Recently, a systematic treatment of correlations in 
highly-coupled systems has been put forward [21], that employs a perturbative scheme (virial 
expansion) in terms of 1/S. This scheme leads to the asymptotic strong-coupling (SC) the- 
ory, which becomes exact in the limit which is complementary to the mean-fleld regime, i.e. 
when S — > oo. The SC theory has been used to study the interaction between planar charged 
walls [21] and shows quantitative agreement with Monte-Carlo simulations for moderate to 
large coupling parameters [22]. The SC mechanism of like-charge attraction qualitatively 
agrees with the structural-correlation scenario [17-20] for large macroion charges [21-23,27]. 
In this letter, we study the effective interaction in a system of two like-charged rods using 
both Molecular Dynamics (MD) simulation methods and the SC theory. The nature of elec- 
trostatic correlations in such a system has been studied in our previous MD simulations [10], 
and exhibits a competition between electrostatic and excluded-volume interactions. Here, we 
focus on the strong-coupling characteristics of attraction between like-charged rods. In the 
following, these numerical results will be compared with the predictions of the SC theory [27]. 

Strong- coupling theory. - In the limit of large coupling parameter S ;3> 1, the canonical 
free energy of a system of fixed macroions with their neutralizing counterions, admits a 
large-coupling expansion as T = J-ifB. + + . . . [27], where the coefficients J^i, J^2, • ■ • 

are expressed in terms of weighted integrals of Mayer functions over counterionic degrees 
of freedom, which are convergent for systems of counterions at charged macroscopic objects 
[21]. The first term in the above expression generates the leading non- vanishing contribution 
to statistical observables such as effective interaction between macroions for S ^ 1; hence 
-^SC = -^i/S is referred to as the SC free energy of the system and as the rescaled 
SC free energy. The SC free energy takes a very simple form in terms of the one-particle 
interaction energies between counterions and macroions (see, e.g., eq. below), and the 
mutual interaction between counterions themselves enters only in the sub-leading corrections 
[21,27]. Physically, this reflects the fact that for S ^ 00, counterions at charged objects are 
surrounded by an extremely large correlation hole and thus, the single-particle term becomes 
the dominant contribution [21-23,27] (see also the discussion). Let us consider two infinitely 
long similar rods of length H and radius R that are located parallel at axial separation D in 
a rectangular box of lateral extension L containing also neutralizing counterions of valency 
q (fig. ^). The rods are impenetrable to counterions and no additional salt is present. It 
is convenient to establish a dimcnsionlcss formalism for such systems by rescaling all length 
scales by the Gouy-Chapman length /i — l/(27r(7^Bfs), such as x = x/fi, etc. [21], where 
(Ts = r/(27ri?) stands for the surface charge density of the rods. Clearly, the rod radius in 
rescaled units, R = coincides with Manning parameter, i.e. R = ^ = qt-eT. (We may 

occasionally use ^ or i? to denote the Manning parameter.) The rescaled rod length, H , is 
related to the total number of counterions, TV, via the global electroneutrality condition as 
H = NE/^. The rescaled SC free energy of this system per rescaled unit length follows as [27] 
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Fig. 1 - a) We consider two identical rods of radius R at axial distance D in an outer box of size 
L, which is rectangular (with square cross-section) in the analytical model and cylindrical in the 
simulations. i?o and a are the Lennard- Jones potential parameters defined in the simulation model, 
b) The rescaled equilibrium surface-to-surface distance of rods, A, = A»/ fi — (D* — 2R)/ fi, obtained 
from the SC theory as a function of Manning parameter, ^ — qisT, for various rescaled box sizes, 
L//i, indicated on the graph. In the infinite-box limit, a continuous unbinding occurs at ~ 2/3. 



H 



dxdy 



(1) 



The integral runs over the volume V available for counterions within the confining box, ex- 
cluding the rods. We defined fi,2 — [{x ± i)/2)^ as the radial distances from the rods 
axes. The first term in eq. is nothing but the bare electrostatic repulsion between the rods. 
The second term involves the leading (energetic and entropic) contributions from counterions, 
which leads to a counterion-mediated attraction between sufficiently highly charged rods [29] . 
This term also reflects the de-condensation process of counterions at low Manning parame- 
ters: For ^ < 1/2, the counterionic integral in eq. 1^ diverges in the infinite- volume limit 
L — > oo, hence the distribution of counterions around the rods and the counterion-mediated 
force vanish, i.e. the rods purely repel each other [25]. For large Manning parameter ^ 3> 1, 
on the other hand, the rescaled SC free energy, J^i, shows a long-ranged attraction and a 
pronounced global minimum at a small axial separation _D* w 2i?, which is nothing but the 
equilibrium axial separation of the rods [27] . For ^ » 1 , the counterionic distribution becomes 
strongly localized in a narrow region between the rods, indicating that the SC attraction is 
accompanied by strong accumulation of counterions in the intervening region. This allows for 
a saddle-point calculation of the counterionic integral in eq. Q , giving the approximate form 
of jFi in the vicinity of its local minimum as [27] 

Ti/H « 6^^ In - 2^ ln(L> - 2R). (2) 

The equilibrium surface-to-surface distance of the rods for large ^ is obtained by minimizing 
expression as A* = — 27? sa 2/3-1-0(1/^), which corresponds, in real units, to a surface 
separation of the order of the Gouy-Chapman length, fi [31]. Figure^ shows the global 
behavior of the rescaled equilibrium surface-to-surface distance, A,, obtained by numerical 
minimization of the full SC free energy, eq. (^3), as a function of f for several box sizes. 
For large ^, attraction is dominant, the rods form a closely-packed bound state, and the 
volume of the bounding box is irrelevant due to strong condensation of counterions. However, 
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for decreasing ^, attraction is weakened and rods eventually unbind in a continuous fashion 
when L oo. For two unconfined rods (solid curve in fig. ^p), the unbinding transition 
occurs at = 2/3 and exhibits an asymptotic power-law behavior as a function of the 
reduced Manning parameter as A* ~ — 2/3) where a w 3/2 [27]. Note that the 
predicted onset of rod- rod attraction, — 2/3, is somewhat larger than the counterion- 
condensation threshold for the combined system of two rods, ^ = 1/2 [26], and smaller than 
the condensation threshold for a single rod, ^ = 1 [24]. This reflects the subtle crossover from 
the case of two coupled rods (when the distance between them is small) to the case of two 
decoupled rods as the axial distance diverges when the unbinding threshold is approached. The 
predicted attractive force between rods, i^rods = —{kBT)dJ-'g,c/dD (in actual units), is found 
from eq. ^ to be inversely proportional to the axial distance, D. This force increases and 
tends to a temperature- independent value upon increasing ^ (or decreasing temperature) [3] , 
which is nothing but the energetic attraction mediated by counterions sandwiched between 
closely-packed rods, i.e. i^rods = — 3e-^T^/(27reeo-D) [27]. The SC attraction for large ^, 
thus, originates from the energetic contributions induced by the ground-state configuration of 
counterions and in this respect, agrees qualitatively with the low-temperature picture for like- 
charge attraction [17-20]. Though for the specific case of two charged rods, the quantitative 
predictions of the SC theory for the (energetic) inter-rod force at large ^, and for the attraction 
threshold differ from the results of the model studied in ref . [20] , which incorporates a discrete 
charge pattern for the rods and obtains attraction for ^ > 2 [20] . 

Simulation method. - To study the interaction between two like-charged rods numeri- 
cally, we have performed extensive Molecular Dynamics simulations making use of a Langevin 
thermostat to drive the system into the canonical state. The geometry of the simulated sys- 
tem is similar to what we have sketched in fig. QJi; the outer confining box is here chosen to 
be cylindrical. Other parameters are defined in the same way unless explicitly mentioned. 
Charged particles interact via bare Coulombic interaction and excluded-volume interactions. 
We apply periodic boundary conditions in the direction parallel to the rods axes. This leads to 
summation of infinite series for Coulomb interaction over all periodic images, which is handled 
numerically using the MMMID summation method [30]. As we are only interested in exam- 
ining electrostatic aspects of the effective rod-rod interaction, excluded-volume interactions 
are considered only between rods and counterions and not between counterions themselves. 
(Extended numerical results with excluded-volume interaction between counterions will be 
presented elsewhere [28].) We employ a shifted Lennard- Jones potential as 



where Rq is an offset used to control the radius of the rods, a defines our basic length scale 
in the simulations, and e « Ik^T. Therefore, R ~ Rq + a may be regarded as the effective 
rod radius and is directly compared with the SC result (with hard-core rods). This is jus- 
tified since in the simulations only a negligible fraction of counterions penetrates the above 
potential to reach radial distances smaller than R. In the simulations, the diameter of the 
(cylindrical) outer box is chosen as 8D. For the final comparison (fig.|3), the theoretical curves 
are also calculated using a similar constraint, though for simplicity and as explained before, 
calculations have been done for a square box of edge size L — 8D. As we have explicitly 
checked, the results are insensitive to small changes in the box size for the considered range of 
Manning parameter (see also fig. QJd) , and the SC predictions can thus be compared with the 
simulations. To obtain the location of the zero-force (equilibrium) point in the force curves. 




otherwise 




a 



(3) 
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Fig. 2 - The equilibrium surface-to-surface distance of two rods in rescaled units, A* = A*//i = 
(-D* — 2R)/iJ., as a function of Manning parameter, ^ — qIbt. Symbols show simulation results for 
(from top): 7rb = 3 (crosses), 9.1 (open circles), 10 (open triangle-downs), 15 (open triangle-ups), 
30 (stars), 40 (open squares), 50 (open diamonds), 60 (filled diamonds). The solid curve is the SC 
prediction obtained by numerical minimization of eq. Q. Dashed curves are guides to the eye. 

we have employed a bisection algorithm followed by a linear regression. The error bars are 
determined from the error of the linear regression using the error propagation method. 

Results and discussion. - To clearly separate the attraction and repulsion regimes of two 
like-charged rods numerically, we study the rescaled equilibrium surface-to-surface separation 
of rods. A*, as a function of Manning parameter, ^, in our simulations. The results will then 
be compared with the SC prediction (fig.n})). To this end, we fix the actual surface-to-surface 
distance between the rods, A = D — 2R, and their linear charge density, r, together with 
the counterion valency, q, but vary the Bjerrum length, ^b, and the actual rod radius, R. By 
doing so, the Gouy-Chapman length associated with this system, fj, — R/{£b<1t), is varied, 
which then allows for determining the equilibrium surface-to-surface separation in rescaled 
units, A* = A^/fi. The results will also depend upon the electrostatic coupling parameter, 
5 — 2TTq^i'^as, which is finite in the simulations and may be written as 

S = eA7RB, (4) 

where 7r,b = '//(''"A) is a dimensionless parameter referred to here as the Rouzina-Bloomfield 
parameter (see the discussion below). The simulation results for the equilibrium surface- 
to-surface separation of the rods are shown in fig. |21 for different 7r,b ranging from 3 to 
60. The data have been obtained using various combinations of fixed ra — 1.0,0.33,0.10, 
q = 1,3, 5, 10 and A/a — 0.5, 1.0, 2.0, 2.5. For a given set of system parameters (r, q and A) 
chosen in the simulations, 7rb is fixed and identifies a single curve. As seen, upon increasing 
the Rouzina-Bloomfield parameter, 7r,b, the equilibrium separation decreases indicative of a 
stronger attraction operating between the rods, and at the same time, the agreement between 
simulation data and the SC prediction (solid curve in fig.|21) becomes progressively better. The 
agreement is quantitative for large 7rb, *-e. 7rb = 50 and 60. (The data set with 7rb = 60 
is obtained using q = 3, ra = 0.1 and A/a = 0.5, and data with 7rb = 50 are obtained using 
5 = 5, TO = 0.1 and A/a — 1.0.) The observed trend for increasing 7rb is associated with the 
increase of the coupling parameter, 2, in the simulated system (see eq. (0))), which eventually 
exhibits the strong-coupling regime. (For instance, for a moderate Manning parameter of 
^ = 3.0, S increases from 18 for cross symbols up to about 100 for filled diamonds.) 
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It is noteworthy that 7rb, as defined above, can be expressed as 7rb = <5/A, where 
5 = q/r is the projected distance between countcrions along a single rod as follows from the 
local electroneutrality condition. Note that 6 roughly gives the correlation hole size around 
counterions at charged surfaces [6,7,10,22,23]. The appearance of strong attractive forces 
for 7RB > 1 was first addressed by Rouzina and Bloomfield [17] for a system of two planar 
charged walls. Analytical and numerical results for planar walls [21-23] have in fact shown 
that the relative magnitude of higher-corrections to the asymptotic SC theory decreases with 
7rb5 and the SC regime at finite S is characterized by the Rouzina-Bloomfield criterion 
7rb > 1- Physically, this corresponds to a large correlation hole size around counterions at 
macroion surfaces as compared with the macroions surface separation, i.e. 6 > A, leading to 
a dominant contribution from single (isolated) counterions [21-23], which is formally obtained 
within the SC scheme [21]. It becomes exceedingly difiicult to perform systematic analytic 
calculations to examine higher-order corrections to SC predictions given the geometry of the 
two-rod system. Nonetheless, the present numerical results clearly indicate a qualitatively 
similar trend for increasing 7r,b in this system and that, the SC regime is characterized by 
Rouzina-Bloomfield criterion. Due to convergence limitations, our numerical investigation so 
far has been limited to the range of ^ > 0.8 and thus do not cover the close vicinity of the 
unbinding transition. (It becomes more difiicult to obtain good data as the rods equilibrium 
distance rapidly increases for small ^.) However, the excellent convergence of the simulation 
results to the asymptotic SC prediction suggests an attraction threshold of about = 2/3 
for two rods in the strong-coupling limit as obtained using the SC theory. Note that for small 
Manning parameters where counterions de-condense (^ < 1/2), electrostatic correlations are 
suppressed due to entropic dilution of the counterionic cloud. Physically, this regime does not 
exhibit strong energetic coupling and higher-order corrections to the SC theory may become 
important. Nevertheless as we showed, the SC theory captures the de-condensation process 
and leads to a qualitatively consistent picture for the whole range of Manning parameters [25]. 
(Note that the predicted onset of de-condensation at ^ = 1/2 and also the bare macroionic 
repulsion for small ^ quantitatively agree with previous findings [26].) How and whether the 
predicted attraction threshold (^c = 2/3) and the unbinding behavior is influenced by the 
de-condensation process of counterions at small ^ remains to be clarified. 

In summary, we have investigated the attraction and repulsion regimes of two like-charged 
rods in terms of Manning parameter, ^, and the coupling parameter, S, by combining nu- 
merical and analytical approaches. For large ^, the rods form a closely-packed bound state 
of small surface-to-surface separation of about the Gouy-Chapman length, fi [31]. The at- 
traction is weakened and the rods unbind for decreasing ^. The attraction threshold is found 
at = 2/3, and the unbinding of rods proceeds in a continuous fashion characterized by a 
power-law for the rods surface separation. A* ~ (^ — 2/3)~", with an exponent a « 3/2. The 
predicted attraction regime for large ^ is accessible experimentally using, e.g., DNA molecules 
(re « 6e/nm) in the presence of multivalent counterions such as spermidine (g = 3), which 
yields ^ « 12 and S w 80. For these values, our results predict attraction with an equilibrium 
surface separation A* ^ /i, where for the DNA-spermidine system /i ^ lA [31]. It should be 
noted, however, that for DNA-like molecules (large r), other factors such as excluded- volume 
interaction between counterions [10,28] and the helical structure of DNA [18, 19] may become 
important and contribute additional components to the total effective force. The predictcid 
continuous unbinding occurs at quite small values of Manning parameter. Experimentally, 
this phenomenon could be studied with weakly charged stiff polymers, such as poly-phenylene 
with a suitably chosen small density of charged side-chains. The effects of salt can be quali- 
tatively accounted for by associating the bounding box size with the screening length: Thus 
only close to the unbinding do we expect addition of salt to matter. The eff'ects of finite chain 
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stiffness and the bundling of many chains constitute interesting apphcations for the future. 
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